################################################################################
##Replication files for
##Moral Individualism in Modern Politics: A New Measure Inspired by Political Theory
##Christopher F. Karpowitz and Kelly D. Patterson
##Perspectives on Politics
################################################################################

library(pacman)
p_load(varhandle, tidyverse, foreign, arm, stargazer, reshape2, data.table, 
       ggplot2, psych, GPArotation, sp, readstata13, plyr, plm, pglm, nFactors,
       margins, rsvg, ggpubr, labelled, weights, jtools, estimatr, xtable, lme4,
       marginaleffects, janitor, RStata, here)

#load original data
##Set working directory to location of replication folder
i_am("KarpowitzPatterson_Perspectives_Replication.R")


###########################################################
##Load datasets
###########################################################
##Setup RStata
options("RStata.StataPath" = "/Applications/Stata/StataSE.app/Contents/MacOS/stata-se")
options("RStata.StataVersion" = 17)

yougov <-read.dta13("YouGov_replication.dta")

covid <- read.dta13("YouGov_Coronavirus_replication.dta")

western <- read.dta13("WesternStates2020_replication.dta")

unc <- read.dta13("COVID UNC Waves 6 and 7_replication.dta")

##GSS Data can be downloaded from this website: https://gss.norc.org/get-the-data/stata
## This is the version we used: GSS 2016-2020 Panel (Release 1a, April 2022)
gss <- read.dta13("gss2020panel_r1a.dta")


###########################################################
##Summary statistics for each study
###########################################################

##Appendix Table A1
stargazer(yougov[c("indiv_self_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "authoritar", "ideo5_new", "pid7", "age", "female", "race_nonwhite", "college_grad", "faminc_quintile")],
          type="latex", out="Tables/yougov_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Partisanship", "Age", "Female", "Nonwhite", "College Graduate", "Family Income (Quintile)") 
)

##Appendix Table A2
covid$pid7_summary <- as.numeric(covid$pid7)
stargazer(covid[c("indiv_moral_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "authoritar", "ideo5_new", "pid7_summary", "age", "female", "race_nonewhite", "college_grad", "faminc_quintile")],
          type="latex", out="Tables/covid_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Partisanship", "Age", "Female", "Nonwhite", "College Graduate", "Family Income (Quintile)") 
)
covid <- subset(covid, select = -c(pid7_summary))

##Appendix Table A3
stargazer(western[c("indiv_moral_irt_group", "indiv_econ", "authoritar", "ideo5_new", "pid7", "age", "female", "race_nonewhite", "college_grad", "faminc_quintile")],
          type="latex", out="Tables/western_summary.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Authoritarianism", "Ideology", "Partisanship", "Age", "Female", "Nonwhite", "College Graduate", "Family Income (Quintile)") 
)

##Appendix Table A4
stargazer(unc[c("indiv_moral_irt_group_w6", "ideo5_new_w6", "pid7_w6", "age_w6", "female_w6", "race_nonwhite_w6", "college_grad_w6", "faminc_quintile_w6")],
          type="latex", out="Tables/unc_summary_w6.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Ideology", "Partisanship", "Age", "Female", "Nonwhite", "College Graduate", "Family Income (Quintile)") 
)

##Appendix Table A5
stargazer(unc[c("indiv_moral_irt_group_w7", "authoritar", "ideo5_new_w7", "pid7_w7", "age_w7", "female_w7", "race_nonwhite_w7", "college_grad_w7", "faminc_quintile_w7")],
          type="latex", out="Tables/unc_summary_w7.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Ideology", "Partisanship", "Age", "Female", "Nonwhite", "College Graduate", "Family Income (Quintile)") 
)


###########################################################
##Variable recoding for analysis
###########################################################
yougov$pid7 <- factor(yougov$pid7,
                      levels = c(1, 2, 3, 4, 5, 6, 7),
                      labels = c("Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Independent", "Lean Republican", "Not very strong Republican", "Strong Republican"))

yougov <- within(yougov, pid7 <- relevel(factor(pid7), ref = 4))

covid <- covid %>% mutate(pid7 = case_when(pid7=="Strong Democrat" ~ 1,
                                           pid7=="Not very strong Democrat" ~ 2,
                                           pid7=="Lean Democrat" ~ 3,
                                           pid7=="Independent" ~ 4,
                                           pid7=="Lean Republican" ~ 5,
                                           pid7=="Not very strong Republican" ~ 6,
                                           pid7=="Strong Republican" ~ 7))
covid$pid7 <- factor(covid$pid7,
                     levels = c(1, 2, 3, 4, 5, 6, 7),
                     labels = c("Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Independent", "Lean Republican", "Not very strong Republican", "Strong Republican"))

covid <- within(covid, pid7 <- relevel(factor(pid7), ref = "Independent"))

western$pid7 <- factor(western$pid7,
                     levels = c(1, 2, 3, 4, 5, 6, 7),
                     labels = c("Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Independent", "Lean Republican", "Not very strong Republican", "Strong Republican"))

western <- within(western, pid7 <- relevel(factor(pid7), ref = 4))


unc$pid7_w6 <- factor(unc$pid7_w6,
                     levels = c(1, 2, 3, 4, 5, 6, 7),
                     labels = c("Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Independent", "Lean Republican", "Not very strong Republican", "Strong Republican"))

unc <- within(unc, pid7_w6 <- relevel(factor(pid7_w6), ref = 4))


unc$pid7_w7 <- factor(unc$pid7_w7,
                      levels = c(1, 2, 3, 4, 5, 6, 7),
                      labels = c("Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Independent", "Lean Republican", "Not very strong Republican", "Strong Republican"))

unc <- within(unc, pid7_w7 <- relevel(factor(pid7_w7), ref = 4))


###########################################################
##Analysis
###########################################################

##Table 1: Authority Referents: Most Important Influence in Helping Make Decisions
##Summary stats of referents
yougov_ref <- rd(100*wpct(yougov$CSED100_2, yougov$weight), digits=1)
covid_ref <- rd(100*wpct(covid$referent, covid$weight), digits=1)
western_ref <- rd(100*wpct(western$referent, western$weight), digits=1)
unc_w6_ref <- rd(100*wpct(unc$referent_w6), digits=1)
unc_w7_ref <- rd(100*wpct(unc$referent_w7), digits=1)

yougov_refcount <- sum(with(yougov, !is.na(CSED100_2)))
covid_refcount <- sum(with(covid, !is.na(referent)))
western_refcount <- sum(with(western, !is.na(referent)))
unc_w6_refcount <- sum(with(unc, !is.na(referent_w6)))
unc_w7_refcount <- sum(with(unc, !is.na(referent_w7)))

refcount<-rbind(yougov_refcount, covid_refcount, western_refcount, unc_w6_refcount, unc_w7_refcount)
refcount <- t(refcount)

ref <- cbind(yougov_ref, covid_ref, western_ref, unc_w6_ref, unc_w7_ref)
ref <- rbind(ref, refcount)
rownames(ref) <- c("Family", "Religion", "Science", "A Good Friend", "Other", "N")

stargazer(ref,
          type="latex", out="Tables/summary_referent.tex", float = FALSE, #style="apsr",
          covariate.labels   = c("\\thead{Referent}", "\\thead{YouGov 2018}", "\\thead{YouGov 2020}", "\\thead{WSS 2020}", "\\thead{UNC 2021}", "\\thead{UNC 2022}")#,
          #add.lines=list(c("N", refcount))#,
          #notes = "\\parbox[t]{.7\\textwidth}{Cell entries are OLS regression coefficients; standard errors in parentheses. $^{*}$p$<$0.10; $^{**}$p$<$0.05; $^{***}$p$<$0.01, two-tailed.}",
          #notes.append = FALSE, # TRUE append the significance levels
          #notes.align = 'l' 
)

###########################################################
##Distribution of moral individualism by referent
###########################################################

##Create kdensity plots of key variables
mdat_2018 <- summarise(yougov, loc.mean=mean(indiv_self_irt_group), loc2.mean=mean(indiv_behave), loc3.mean=mean(indiv_econ), loc4.mean=mean(authoritar), loc5.mean=mean(collect_all), loc6.mean=mean(indiv_horiz), loc7.mean=mean(indiv_vert))
yougov$ref <- cut(yougov$CSED100, breaks=c(0, 1, 2, 3, 4, 5, 6, 7, 8), labels=c("Religion", "Family", "Science", "Teacher", "Work Colleague", "Public", "Friend", "Media"))
yougov$ref2[yougov$CSED100==2] <- "Family"
yougov$ref2[yougov$CSED100==1] <- "Religion"
yougov$ref2[yougov$CSED100==3] <- "Science"
yougov$ref2[yougov$CSED100==7] <- "Friend"
yougov$ref2[yougov$CSED100==4] <- "Other"
yougov$ref2[yougov$CSED100==5] <- "Other"
yougov$ref2[yougov$CSED100==6] <- "Other"
yougov$ref2[yougov$CSED100==8] <- "Other"
mdat2_2018 <- ddply(yougov, "ref2", summarise, loc.mean=mean(indiv_self_irt_group), loc3.mean=mean(indiv_econ), loc4.mean=mean(authoritar), loc6.mean=mean(indiv_horiz), loc7.mean=mean(indiv_vert))


##Appendix Figure A1, Panel A
##Distribution of moral individualism by referent -- 2018
moral_indiv_referent_2018 <- ggplot(yougov, aes(x=indiv_self_irt_group, colour = ref2, fill = ref2)) +
  geom_density(alpha=0.6)+
  geom_vline(data=mdat2_2018, aes(xintercept=loc.mean,  colour=ref2),
             linetype="dashed", size=1)+
  scale_fill_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), breaks=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") + 
  scale_color_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), breaks=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") +
  xlab("Moral Individualism") +
  theme_minimal()
moral_indiv_referent_2018
#ggsave("dist_referent_2018.pdf", width = 10.57, height = 9.55, dpi = 320)



##Figure 1: Distribution of moral individualism by referent, YouGov 2020
##Create kdensity plots of key variables
covid$ref2 <- as.factor(covid$referent)
mdat2_2020 <- ddply(covid, "ref2", summarise, loc.mean=mean(indiv_moral_irt_group), loc2.mean=mean(indiv_moral_all), loc3.mean=mean(indiv_econ), loc4.mean=mean(authoritar), loc5.mean=mean(indiv_goertz), loc6.mean=mean(indiv_horiz), loc7.mean=mean(indiv_vert))

##2020 distribution
moral_indiv_referent_2020 <- ggplot(covid, aes(x=indiv_moral_irt_group, colour = ref2, fill = ref2)) +
  geom_density(alpha=0.6)+
  geom_vline(data=mdat2_2020, aes(xintercept=loc.mean,  colour=ref2),
             linetype="dashed", size=1)+
  scale_fill_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") + 
  scale_color_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") +
  xlab("Moral Individualism") +
  theme_minimal()
moral_indiv_referent_2020
#ggsave("Figures/moral_indiv_referent_2020.pdf", width = 10.57, height = 9.55, dpi = 320)

##Figure A1, Panel B
##WSS 2020 distribution
western$ref2 <- as.factor(western$referent)
mdat2_western <- ddply(western, "ref2", summarise, loc.mean=mean(indiv_moral_irt_group), loc3.mean=mean(indiv_econ), loc4.mean=mean(authoritar))

moral_indiv_referent_western <- ggplot(western, aes(x=indiv_moral_irt_group, colour = ref2, fill = ref2)) +
  geom_density(alpha=0.6)+
  geom_vline(data=mdat2_western, aes(xintercept=loc.mean,  colour=ref2),
             linetype="dashed", size=1)+
  scale_fill_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") + 
  scale_color_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") +
  xlab("Moral Individualism") +
  theme_minimal()
moral_indiv_referent_western
#ggsave("dist_referent_western.pdf", width = 10.57, height = 9.55, dpi = 320)


##UNC 2021 and UNC 2022 distribution
unc$ref2_w6 <- as.factor(unc$referent_w6)
mdat2_unc_w6 <- ddply(unc, "ref2_w6", summarise, loc.mean=mean(indiv_moral_irt_group_w6))

unc$ref2_w7 <- as.factor(unc$referent_w7)
mdat2_unc_w7 <- ddply(unc, "ref2_w7", summarise, loc.mean=mean(indiv_moral_irt_group_w7))

moral_indiv_referent_unc21 <- ggplot(subset(unc, !(is.na(referent_w6))), aes(x=indiv_moral_irt_group_w6, colour = ref2_w6, fill = ref2_w6)) +
  geom_density(alpha=0.6)+
  geom_vline(data=mdat2_unc_w6, aes(xintercept=loc.mean,  colour=ref2_w6),
             linetype="dashed", size=1)+
  scale_fill_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") + 
  scale_color_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") +
  xlab("Moral Individualism") +
  theme_minimal()
moral_indiv_referent_unc21

moral_indiv_referent_unc22 <- ggplot(unc, aes(x=indiv_moral_irt_group_w7, colour = ref2_w7, fill = ref2_w7)) +
  geom_density(alpha=0.6)+
  geom_vline(data=mdat2_unc_w7, aes(xintercept=loc.mean,  colour=ref2_w7),
             linetype="dashed", size=1)+
  scale_fill_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") + 
  scale_color_manual(values=c("#d53e4f", "#fc8d59", "#fee08b", "#99d594", "#3288bd"), labels=c("Family", "Religion", "Science", "Friend", "Other"), name="Authority") +
  xlab("Moral Individualism") +
  theme_minimal()
moral_indiv_referent_unc22

##Figure A1: Moral Individualism Distributions by Referent
##Combine desnsity plots
ggarrange(moral_indiv_referent_2018, moral_indiv_referent_western, moral_indiv_referent_unc21, moral_indiv_referent_unc22,
          ncol =2, nrow = 2,
          labels = "AUTO",
          #labels = c("2018 YouGov", "2020 Western States", "2021 UNC", "2022 UNC"),
          hjust = 0,
          vjust = 1,
          common.legend = TRUE,
          legend = "bottom")


###Factor Analysis
#New labels -- 2018 National
yougov$MORAL1 <- yougov$CSED102
yougov$MORAL2 <- yougov$CSED104
yougov$MORAL3 <- yougov$CSED105
yougov$MORAL5 <- yougov$CSED106
yougov$MORAL6 <- yougov$CSED107
yougov$MORAL7 <- yougov$CSED108
yougov$MORAL8 <- yougov$CSED127
yougov$MORAL9 <- yougov$CSED128
yougov$MORAL10 <- yougov$CSED130

yougov$ECON1 <- yougov$CSED200
yougov$ECON2 <- yougov$CSED201
yougov$ECON3 <- yougov$CSED202
yougov$ECON4 <- yougov$CSED203
yougov$ECON5 <- yougov$CSED204
yougov$ECON6 <- yougov$CSED205

yougov$HORIZ1 <- yougov$CSED251
yougov$HORIZ2 <- yougov$CSED252
yougov$HORIZ3 <- yougov$CSED253
yougov$HORIZ4 <- yougov$CSED254

yougov$VERT1 <- yougov$CSED255
yougov$VERT2 <- yougov$CSED256
yougov$VERT3 <- yougov$CSED257
yougov$VERT4 <- yougov$CSED258

yougov$AUTH1 <- yougov$CSED400A
yougov$AUTH2 <- yougov$CSED400B
yougov$AUTH3 <- yougov$CSED400C
yougov$AUTH4 <- yougov$CSED400D

factorvars_2018 <- c("MORAL1", "MORAL2", "MORAL3", "MORAL5", "MORAL6", "MORAL7", "MORAL8", "MORAL9", "MORAL10", "ECON1", "ECON2", "ECON3", "ECON4", "ECON5", "ECON6", "HORIZ1", "HORIZ2", "HORIZ3", "HORIZ4", "VERT1", "VERT2", "VERT3", "VERT4", "AUTH1", "AUTH2", "AUTH3", "AUTH4")
yougov_factor <-yougov[factorvars_2018]

##2018
##Appendix Figure A1, Panel B
##Parallel analysis
parallel_2018 <-fa.parallel(yougov_factor, fm='minres', fa='fa')

##Appendix Figure A3: Factor Loadings, 2018 YouGov Survey
fivefactor_2018 <-fa(yougov_factor, nfactors=5, rotate="oblimin", fm="minres")
fivefactor_2018
print(fivefactor_2018$loadings, digits=2, cutoff=.2, sort=TRUE)
fa.diagram(fivefactor_2018, cut=0.25, main="2018 Factor Loadings", size=1)


#New labels -- 2020 National
covid$MORAL1 <- covid$Q4
covid$MORAL2 <- covid$Q5
covid$MORAL3 <- covid$Q6
covid$MORAL4 <- covid$Q7
covid$MORAL5 <- covid$Q8
covid$MORAL6 <- covid$Q9
covid$MORAL7 <- covid$Q10
covid$MORAL8 <- covid$Q11
covid$MORAL9 <- covid$Q12
covid$MORAL10 <- covid$Q13

covid$ECON1 <- covid$Q14_1
covid$ECON2 <- covid$Q14_2
covid$ECON3 <- covid$Q14_3
covid$ECON4 <- covid$Q14_4
covid$ECON5 <- covid$Q14_5
covid$ECON6 <- covid$Q14_6

covid$HORIZ1 <- covid$Q20_1
covid$HORIZ2 <- covid$Q20_2
covid$HORIZ3 <- covid$Q20_3
covid$HORIZ4 <- covid$Q20_4

covid$VERT1 <- covid$Q21_1
covid$VERT2 <- covid$Q21_2
covid$VERT3 <- covid$Q21_3
covid$VERT4 <- covid$Q21_4

covid$AUTH1 <- covid$Q16
covid$AUTH2 <- covid$Q17
covid$AUTH3 <- covid$Q18
covid$AUTH4 <- covid$Q19


##Factor analysis without obedience measures
factorvars <- c("MORAL1", "MORAL2", "MORAL3", "MORAL4", "MORAL5", "MORAL6", "MORAL7", "MORAL8", "MORAL9", "MORAL10", "ECON1", "ECON2", "ECON3", "ECON4", "ECON5", "ECON6", "HORIZ1", "HORIZ2", "HORIZ3", "HORIZ4", "VERT1", "VERT2", "VERT3", "VERT4", "AUTH1", "AUTH2", "AUTH3", "AUTH4")

covid_factor <-covid[factorvars]

covid_factor <- covid_factor %>%
  mutate_if(is.factor, as.numeric)



##2020
##Parallel analysis
##Figure A2, Panel A 
parallel <-fa.parallel(covid_factor, fm='minres', fa='fa')

##Figure 2: Factor Analysis of Key Measures, YouGov 2020
fivefactor <-fa(covid_factor, nfactors=5, rotate="oblimin", fm="minres")
print(fivefactor)
print(fivefactor$loadings, digits=2, cutoff=.2, sort=TRUE)
fa.diagram(fivefactor, cut=0.25, main="Factor Loadings")
#biplot(fivefactor)
fa.parallel(covid_factor)

##Evidence for patterns with additional factors
sixfactor <-fa(covid_factor, nfactors=6, rotate="oblimin", fm="minres")
print(sixfactor$loadings, digits=2, cutoff=.2, sort=TRUE)

sevenfactor <-fa(covid_factor, nfactors=7, rotate="oblimin", fm="minres")
print(sevenfactor$loadings, digits=2, cutoff=.2, sort=TRUE)


##Figure A4: Correlations between Key Measures, 2018 YouGov Survey
##Correlation matrix
corrs_2018 <- c("indiv_self_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "authoritar", "ideo5_new")
corrs <- c("indiv_moral_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "authoritar", "ideo5_new")
corrs_western <- c("indiv_moral_irt_group", "indiv_econ", "authoritar", "ideo5_new")

corr_data_2018 <- yougov[corrs_2018]
lowerCor(corr_data_2018)
pairs.panels(corr_data_2018,pch='.', cex.labels=0.6, labels=c("Moral", "Economic", "Horizontal", "Vertical", "Authoritarianism", "Ideology"))


##Figure 3: Correlations between Key Measures
##YouGov 2020
corr_data_2020 <- covid[corrs]
lowerCor(corr_data_2020)
#plot(threefactor)
pairs.panels(corr_data_2020,pch='.', cex.labels=0.6, labels=c("Moral", "Economic", "Horizontal", "Vertical", "Authoritarianism", "Ideology"))

corr_data_western <- western[corrs_western]
lowerCor(corr_data_western)
#plot(threefactor)
pairs.panels(corr_data_western,pch='.', cex.labels=0.6, labels=c("Moral", "Economic", "Authoritarianism", "Ideology"))


##Table A8: Demographic Characteristis of Individualists, YouGov 2018
##YouGov 2018
demog_moral_yougov <- lm(indiv_self_irt_group ~ age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights=weight)
summ(demog_moral_yougov)
demog_econ_yougov <- lm(indiv_econ ~ age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights=weight)
summ(demog_econ_yougov)
demog_horiz_yougov <- lm(indiv_horiz ~ age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights=weight)
summ(demog_horiz_yougov)
demog_vert_yougov <- lm(indiv_vert ~ age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights=weight)
summ(demog_vert_yougov)

stargazer(demog_moral_yougov, demog_econ_yougov, demog_horiz_yougov, demog_vert_yougov, 
          type="latex", out="Tables/demog_yougov.tex", float = FALSE, #style="apsr",
          column.labels   = c("\\thead{Moral}", "\\thead{Economic}", "\\thead{Horizontal}", "\\thead{Vertical}"),
          #column.separate = c(2, 2), 
          covariate.labels = c("Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq"), dep.var.labels=c("Determinants of Individualism"), no.space = FALSE,
          omit.table.layout = "n",
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "Yes"))
)


##Table A7: Demographic Characteristis of INdividualists, YouGov 2020
demog_moral_covid <- lm(indiv_moral_irt_group ~ age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights=weight)
summ(demog_moral_covid)
demog_econ_covid <- lm(indiv_econ ~ age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights=weight)
summ(demog_econ_covid)
demog_horiz_covid <- lm(indiv_horiz ~ age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights=weight)
summ(demog_horiz_covid)
demog_vert_covid <- lm(indiv_vert ~ age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights=weight)
summ(demog_vert_covid)

stargazer(demog_moral_covid, demog_econ_covid, demog_horiz_covid, demog_vert_covid, 
          type="latex", out="Tables/demog_covid.tex", float = FALSE, #style="apsr",
          column.labels   = c("\\thead{Moral}", "\\thead{Economic}", "\\thead{Horizontal}", "\\thead{Vertical}"),
          #column.separate = c(2, 2), 
          covariate.labels = c("Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq"), dep.var.labels=c("Determinants of Individualism"), no.space = FALSE,
          omit.table.layout = "n",
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "Yes"))
)


##Western
demog_moral_western <- lm(indiv_moral_irt_group ~ age + female + race_nonewhite + college_grad + faminc_quintile, data=western, weights=weight)
summ(demog_moral_western)
demog_econ_western <- lm(indiv_econ ~ age + female + race_nonewhite + college_grad + faminc_quintile, data=western, weights=weight)
summ(demog_econ_western)

##UNC
var_label(unc) <- list(age_w6 = "Age",
                   female_w6 = "Female",
                   race_nonwhite_w6 = "Nonwhite",
                   college_grad_w6 = "College Graduate",
                   faminc_quintile_w6 = "Family Income")

demog_moral_unc_w6 <- lm(indiv_moral_irt_group_w6 ~ as.numeric(age_w6) + female_w6 + race_nonwhite_w6 + college_grad_w6 + faminc_quintile_w6, data=unc)
summ(demog_moral_unc_w6)
demog_moral_unc_w7 <- lm(indiv_moral_irt_group_w7 ~ as.numeric(age_w7) + female_w7 + race_nonwhite_w7 + college_grad_w7 + faminc_quintile_w7, data=unc)
summ(demog_moral_unc_w7)

##Table A9: Demographic Characteristics of Individualists, Western States 2020 and UnC 2021 and 2022 Surveys
##Table fromatting updated by hand, but all coefficients are correct
stargazer(demog_moral_western, demog_econ_western, demog_moral_unc_w6, demog_moral_unc_w7, 
          type="latex", out="Tables/demog_covid_western_unc.tex", float = FALSE, #style="apsr",
          column.labels   = c("\\thead{Moral \\ Western 2020}", "\\thead{Economic \\ Western 2020}", "\\thead{Moral \\ UNC 2021}", "\\thead{Moral \\ UNC 2022}"),
          #column.separate = c(2, 2), 
          covariate.labels = c("Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq"), dep.var.labels=c("Determinants of Individualism"), no.space = FALSE,
          omit.table.layout = "n",
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "Yes"))
)


##Stability over time
##Table A6: Stability of Referent Choices, UNC 2021 and UNC 2022
##Referent stability
ref_names <- c("Family", "Religion", "Science", "Friend", "Other")
unc$ref_value_w6 <- factor(unc$referent_w6,
                           levels = c(1,2,3,4,5),
                           labels = ref_names)
unc$ref_value_w7 <-  factor(unc$referent_w7,
                            levels = c(1,2,3,4,5),
                            labels = ref_names)

ref_stable <- round(prop.table(table(unc$ref_value_w6, unc$ref_value_w7), 1)*100,2)

print(xtable(ref_stable), booktabs=T, file="Tables/ref_stable.tex")

unc$referent_change <- ifelse(unc$referent_w6!=unc$referent_w7, 1, 0)

table(unc$ref_value_w6[!is.na(unc$referent_change)])

moral_all_corr <- round(cor(unc$indiv_moral_irt_group_w6, unc$indiv_moral_irt_group_w7,  method = "pearson", use="pairwise.complete.obs"), 2)
moral_sameref_corr <- round(cor(unc$indiv_moral_irt_group_w6[unc$referent_change==0], unc$indiv_moral_irt_group_w7[unc$referent_change==0],  method = "pearson", use="pairwise.complete.obs"), 2)
moral_difref_corr <- round(cor(unc$indiv_moral_irt_group_w6[unc$referent_change==1], unc$indiv_moral_irt_group_w7[unc$referent_change==1],  method = "pearson", use="pairwise.complete.obs"), 2)

cor.test(unc$indiv_moral_irt_group_w6, unc$indiv_moral_irt_group_w7,  method = "pearson", use="pairwise.complete.obs")
cor.test(unc$indiv_moral_irt_group_w6[unc$referent_change==0], unc$indiv_moral_irt_group_w7[unc$referent_change==0],  method = "pearson", use="pairwise.complete.obs")
cor.test(unc$indiv_moral_irt_group_w6[unc$referent_change==1], unc$indiv_moral_irt_group_w7[unc$referent_change==1],  method = "pearson", use="pairwise.complete.obs")

unc$ref_change <- ifelse(unc$ref_value_w6!=unc$ref_value_w7, 1, 0)
ref_stable <- round(prop.table(table(unc$ref_change))*100, 0)[1]


##Comparison to GSS
gss$party_change <- ifelse(gss$partyid_1b!=gss$partyid_2, 1, 0)
party_stable <- round(prop.table(table(gss$party_change))*100, 0)[1]

gss$partyid_1b_num <- as.numeric(gss$partyid_1b)
gss$partyid_1b_num[gss$partyid_1b_num==8] <- 4

gss$partyid_2_num <- as.numeric(gss$partyid_2)
gss$partyid_2_num[gss$partyid_2_num==8] <- 4

party_corr <- round(cor(gss$partyid_1b_num, gss$partyid_2_num, method = "pearson", use="pairwise.complete.obs"), 2)

gss <- gss %>% mutate(pid3_1b = case_when(partyid_1b=="strong democrat" ~ 1,
                                          partyid_1b=="not very strong democrat" ~ 1,
                                          partyid_1b=="independent, close to democrat" ~ 1,
                                          partyid_1b=="independent (neither, no response)" ~ 2,
                                          partyid_1b=="independent, close to republican" ~ 3,
                                          partyid_1b=="not very strong republican" ~ 3,
                                          partyid_1b=="strong republican" ~ 3,
                                          partyid_1b=="other party" ~ 2))

gss <- gss %>% mutate(pid3_2 = case_when(partyid_2=="strong democrat" ~ 1,
                                          partyid_2=="not very strong democrat" ~ 1,
                                          partyid_2=="independent, close to democrat" ~ 1,
                                          partyid_2=="independent (neither, no response)" ~ 2,
                                          partyid_2=="independent, close to republican" ~ 3,
                                          partyid_2=="not very strong republican" ~ 3,
                                          partyid_2=="strong republican" ~ 3,
                                          partyid_2=="other party" ~ 2))

gss$party_change2 <- ifelse(gss$pid3_1b!=gss$pid3_2, 1, 0)
party_stable2 <- round(prop.table(table(gss$party_change2))*100, 0)[1]


gss$ideo_change <- ifelse(gss$polviews_1b!=gss$polviews_2, 1, 0)
ideo_stable <- round(prop.table(table(gss$ideo_change))*100, 0)[1]

gss <- gss %>% mutate(ideo3_1b = case_when(polviews_1b=="extremely liberal" ~ 1,
                                           polviews_1b=="liberal" ~ 1,
                                           polviews_1b=="slightly liberal" ~ 1,
                                           polviews_1b=="moderate, middle of the road" ~ 2,
                                           polviews_1b=="slightly conservative" ~ 3,
                                           polviews_1b=="conservative" ~ 3,
                                           polviews_1b=="extremely conservative" ~ 3))

gss <- gss %>% mutate(ideo3_2 = case_when(polviews_2=="extremely liberal" ~ 1,
                                          polviews_2=="liberal" ~ 1,
                                          polviews_2=="slightly liberal" ~ 1,
                                          polviews_2=="moderate, middle of the road" ~ 2,
                                          polviews_2=="slightly conservative" ~ 3,
                                          polviews_2=="conservative" ~ 3,
                                          polviews_2=="extremely conservative" ~ 3))

gss$ideo_change2 <- ifelse(gss$ideo3_1b!=gss$ideo3_2, 1, 0)
ideo_stable2 <- round(prop.table(table(gss$ideo_change2))*100, 0)[1]

ideo_corr <- round(cor(as.numeric(gss$polviews_1b), as.numeric(gss$polviews_2), method = "pearson", use="pairwise.complete.obs"), 2)


##Strength of belief in God
gss$god_change <- ifelse(gss$god_1b!=gss$god_2, 1, 0)
god_stable <- round(prop.table(table(gss$god_change))*100, 0)[1]
god_corr <- round(cor(as.numeric(gss$god_1b), as.numeric(gss$god_2), method = "pearson", use="pairwise.complete.obs"), 2)

##Economic class
gss$class_change <- ifelse(gss$class_1b!=gss$class_2, 1, 0)
class_stable <- round(prop.table(table(gss$class_change))*100, 0)[1]
class_corr <- round(cor(as.numeric(gss$class_1b), as.numeric(gss$class_2), method = "pearson", use="pairwise.complete.obs"), 2)

##Belief in papal infallibility
gss$pope_change <- ifelse(gss$popespks_1b!=gss$popespks_2, 1, 0)
pope_stable <- round(prop.table(table(gss$pope_change))*100, 0)[1]
pope_corr <- round(cor(as.numeric(gss$popespks_1b), as.numeric(gss$popespks_2), method = "pearson", use="pairwise.complete.obs"), 2)

##Racial resentment
gss$resent_change <- ifelse(gss$wrkwayup_1b!=gss$wrkwayup_2, 1, 0)
resent_stable <- round(prop.table(table(gss$resent_change))*100, 0)[1]
resent_corr <- round(cor(as.numeric(gss$wrkwayup_1b), as.numeric(gss$wrkwayup_2), method = "pearson", use="pairwise.complete.obs"), 2)

##Government aid to the poor
gss$help_change <- ifelse(gss$helppoor_1b!=gss$helppoor_2, 1, 0)
help_stable <- round(prop.table(table(gss$help_change))*100, 0)[1]
help_corr <- round(cor(as.numeric(gss$helppoor_1b), as.numeric(gss$helppoor_2), method = "pearson", use="pairwise.complete.obs"), 2)


##Stability
##Percentage stable over time
stable_var <- c("Moral Authority Referent", "Political Party (3 point)", "Ideology (3 point)", "Economic Class (4 point)", "Strength of Belief in God (6 point)", "Views of Papal Infallibility (5 point)", "Racial Resentment (5 point)", "Government Aid to Poor (5 point)")
stable_val <- c(ref_stable, party_stable2, ideo_stable2, class_stable, god_stable, pope_stable, resent_stable, help_stable)

stable <- data.frame(variable = stable_var,
                     value = stable_val)
print(xtable(stable, digits=0), booktabs=T, file="Tables/stable_comp.tex")

##Table 3: Correlation between t1 and t2
##Correlation between t1 and t2
corr_var <- c("All Respondents", "Chose Same Referent", "Chose Different Referent", "Political Party (7 point)", "Ideology (7 point)", "Economic Class", "Strength of Belief in God", "Views of Papal Infallibility", "Racial Resentment", "Government Aid to Poor")
corr_val <- c(moral_all_corr, moral_sameref_corr, moral_difref_corr, party_corr, ideo_corr, class_corr, god_corr, pope_corr, resent_corr, help_corr)
corr_time <- data.frame(variable = corr_var,
                        value = corr_val)
print(xtable(corr_time, digits=2), booktabs=T, file="Tables/corr_time_comp.tex")


##Liberty vs. other values
free1_ols <- lm(liberty_security ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(free1_ols)
free1_random <- lmer(liberty_security ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(free1_random)
##Assign pseudo r-squared value. Need to figure out how to grab this from the summ command.
free1_random_pseudor2 <- 0.25

free2_ols <- lm(freedom_death ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(free2_ols)
free2_random <- lmer(freedom_death ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(free2_random)
free2_random_psuedor2 <- 0.33


##Table 4: Attitudes about Freedom
stargazer(free1_ols, free1_random, free2_ols, free2_random, 
          type="latex", out="Tables/free.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Demographic Controls?", "Yes", "Yes", "Yes", "Yes"),
                          c("Survey Weights", "Yes", "No", "Yes", "No"),
                          c("Pseudo R$^{2}$", "", free1_random_pseudor2, "", free2_random_psuedor2))
)

##Table A10: Attitudes about Freedom: Full Results with Demographic Controls
stargazer(free1_ols, free1_random, free2_ols, free2_random, 
          type="latex", out="Tables/free_full.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", free1_random_pseudor2, "", free2_random_psuedor2))
)


##Forced Choice Questions
##Reopen as soon as possible
reopen1 <- glm(reopen ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
              weights = weight,
              family = binomial(link = "probit"), 
              data = covid)
summ(reopen1)
reopen1_pseudor2 <- 0.46
reopen2 <- glmer(reopen ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent),
                family = binomial(link = "probit"), 
                data = covid)
summ(reopen2)
reopen2_pseudor2 <- 0.50

##Price gouging
price_gouge1 <- glm(price_gouge ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
               weights = weight,
               family = binomial(link = "probit"), 
               data = covid)
summ(price_gouge1)
price_gouge1_pseudor2 <- 0.20
price_gouge2 <- glmer(price_gouge ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent),
                 family = binomial(link = "probit"), 
                 data = covid)
summ(price_gouge2)
price_gouge2_pseudor2 <- 0.32

##Support officials
listen_officials1 <- glm(listen_officials ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                    weights = weight,
                    family = binomial(link = "probit"), 
                    data = covid)
summ(listen_officials1)
listen_officials1_pseudor2 <- 0.30
listen_officials2 <- glmer(listen_officials ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent),
                      family = binomial(link = "probit"), 
                      data = covid)
summ(listen_officials2)
listen_officials2_pseudor2 <- 0.38

##Figure 4: Effects of Individualism on Forced Choice Questions
##Produce marginal effects for key IVs for each model
marg_reopen1 <- summary(marginaleffects(reopen1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_reopen1$type <- "reopen1"
marg_price_gouge1 <- summary(marginaleffects(price_gouge1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_price_gouge1$type <- "price_gouge1"
marg_listen_officials1 <- summary(marginaleffects(listen_officials1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_listen_officials1$type <- "listen_officials1"

marg_forced <- rbind(marg_reopen1, marg_price_gouge1, marg_listen_officials1)
rm(marg_reopen1, marg_price_gouge1, marg_listen_officials1)
marg_forced <- marg_forced %>% mutate(measure = case_when(term=="indiv_moral_irt_group" ~ "Moral",
                                          term=="indiv_econ" ~ "Economic",
                                          term=="indiv_horiz" ~ "Horizontal",
                                          term=="indiv_vert" ~ "Vertical",
                                          term=="ideo5_new" ~ "Ideology",
                                          term=="pid7" ~ "Partisanship"))
marg_forced <- marg_forced %>% mutate(dv = case_when(type=="reopen1" ~ "Reopen Economy",
                                                     type=="price_gouge1" ~ "Price Gouging",
                                                     type=="listen_officials1" ~ "Support Officials"))
marg_forced$measure <- factor(marg_forced$measure, levels=c("Moral", "Economic", "Horizontal", "Vertical", "Ideology", "Partisanship"))

fig_forced_choice <- ggplot(data = marg_forced, aes(x=dv, y=estimate, ymin=estimate-std.error*qnorm(.975), ymax=estimate+std.error*qnorm(.975), shape=as.factor(measure))) + 
  expand_limits(y=c(-.50,.50)) + 
  geom_hline(aes(yintercept=0), linetype="dashed", color="red")+
  geom_point(size=4, stat="identity", position=position_dodge(0.4)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(0.4)) +
  ylab("Change in Probability of Support") +
  xlab("Dependent Variable") +
  coord_flip()+
  facet_wrap(~measure)+
  theme_bw()+
  theme(legend.position = "none")
fig_forced_choice
#ggsave("Figures/fig_forced_choice.pdf", width = 10.57, height = 9.55, dpi = 320)

##Tabla A11: Forced Choice Probit/Random Effects
#Supporting table
stargazer(listen_officials1, listen_officials2, reopen1, reopen2, price_gouge1, price_gouge2, 
          type="latex", out="Tables/forced_full.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Support Officials}", "\\thead{Reopen Economy}", "\\thead{Price Gouging}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", listen_officials1_pseudor2, listen_officials2_pseudor2, reopen1_pseudor2, reopen2_pseudor2, price_gouge1_pseudor2, price_gouge2_pseudor2))
)


##Support for official decisions
##Require people to stay at home
require1 <- lm(as.numeric(q43_3) ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(require1)
require2 <- lmer(as.numeric(q43_3) ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(require2)
require2_pseudor2 <- 0.34

##Close schools
close_schools1 <- lm(as.numeric(q43_5) ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(close_schools1)
close_schools2 <- lmer(as.numeric(q43_5) ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(close_schools2)
close_schools2_pseudor2 <- 0.26

##Close non-essential businesses
close_business1 <- lm(as.numeric(q43_7) ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(close_business1)
close_business2 <- lmer(as.numeric(q43_7) ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(relevel(pid7, ref="Independent")) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(close_business2)
close_business2_pseudor2 <- 0.36


##Table A12: Support for Offiical Decisions (OLS/Random Effects), Full Results with Demographic Controls
stargazer(require1, require2, close_schools1, close_schools2, close_business1, close_business2, 
          type="latex", out="Tables/support_full.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Stay-at-Home \\Orders}", "\\thead{Close \\Schools}", "\\thead{Close \\Businesses}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", require2_pseudor2, "", close_schools2_pseudor2, "", close_business2_pseudor2))
)


##Vignettes
#Hand sanitizer vignette
hand1 <- lm(Q38 ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(hand1)
hand2 <- lmer(Q38 ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(hand2)
hand2_pseudor2 <- 0.08

#Mask donation vignette
mask1 <- lm(Q39 ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(mask1)
mask2 <- lmer(Q39 ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(mask2)
mask2_pseudor2 <- 0.08

#Would not go out of your way to purchase goods to donate
covid$donate <- 6-as.numeric(covid$Q40)

donate1 <- lm(donate ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(donate1)
donate2 <- lmer(donate ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(donate2)
donate2_pseudor2 <- 0.08


##Grab coefficients and standard errors from each model and produce new dataset
hand1_coefs <- data.frame(coef(summary(hand1))[c("indiv_moral_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "ideo5_new", "as.factor(pid7)Strong Republican"), 1:2])
hand1_coefs$dv <- "Hand Sanitizer"
hand1_coefs <- hand1_coefs %>% rownames_to_column(var="term")
mask1_coefs <- data.frame(coef(summary(mask1))[c("indiv_moral_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "ideo5_new", "as.factor(pid7)Strong Republican"), 1:2])
mask1_coefs$dv <- "Donate Masks"
mask1_coefs <- mask1_coefs %>% rownames_to_column(var="term")
donate1_coefs <- data.frame(coef(summary(donate1))[c("indiv_moral_irt_group", "indiv_econ", "indiv_horiz", "indiv_vert", "ideo5_new", "as.factor(pid7)Strong Republican"), 1:2])
donate1_coefs$dv <- "Homeless Donation"
donate1_coefs <- donate1_coefs %>% rownames_to_column(var="term")
vignette_coefs <- rbind(hand1_coefs, mask1_coefs, donate1_coefs)
rm(hand1_coefs, mask1_coefs, donate1_coefs)
vignette_coefs <- vignette_coefs %>%
  clean_names()
vignette_coefs <- vignette_coefs %>% mutate(measure = case_when(term=="indiv_moral_irt_group" ~ "Moral",
                                                          term=="indiv_econ" ~ "Economic",
                                                          term=="indiv_horiz" ~ "Horizontal",
                                                          term=="indiv_vert" ~ "Vertical",
                                                          term=="ideo5_new" ~ "Ideology",
                                                          term=="as.factor(pid7)Strong Republican" ~ "Partisanship"))
vignette_coefs$measure <- factor(vignette_coefs$measure, levels=c("Moral", "Economic", "Horizontal", "Vertical", "Ideology", "Partisanship"))
vignette_coefs$dv <- factor(vignette_coefs$dv, levels=c("Homeless Donation", "Hand Sanitizer", "Donate Masks"))


##Figure 5: Effects of Individualism on Responses to Vignettes
fig_vignette <- ggplot(data = vignette_coefs, aes(x=dv, y=estimate, ymin=estimate-std_error*qnorm(.975), ymax=estimate+std_error*qnorm(.975), shape=as.factor(measure))) + 
  expand_limits(y=c(-.50,.50)) + 
  geom_hline(aes(yintercept=0), linetype="dashed", color="red")+
  geom_point(size=4, stat="identity", position=position_dodge(0.4)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(0.4)) +
  ylab("Effect of Independent Variable") +
  xlab("Dependent Variable") +
  coord_flip()+
  facet_wrap(~measure)+
  theme_bw()+
  theme(legend.position = "none")
fig_vignette
#ggsave("Figures/fig_vignette.pdf", width = 10.57, height = 9.55, dpi = 320)

#Supporting Table
#Table A13: Vignettes (OLS/Random Effects)
stargazer(mask1, mask2, hand1, hand2, donate1, donate2, 
          type="latex", out="Tables/vignettes_full.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Donate Masks}", "\\thead{Hand Sanitizer}", "\\thead{Homeless Donation}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", mask2_pseudor2, "", hand2_pseudor2, "", donate2_pseudor2))
)

##Collective action
covid_neighbors1 <- glm(covid_neighbors ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                       weights = weight,
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_neighbors1)
covid_friends1 <- glm(covid_friends ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                        weights = weight,
                        family = binomial(link = "probit"), 
                        data = covid)
summ(covid_friends1)
covid_family1 <- glm(covid_family ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                        weights = weight,
                        family = binomial(link = "probit"), 
                        data = covid)
summ(covid_family1)
covid_online1 <- glm(covid_online ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                     weights = weight,
                     family = binomial(link = "probit"), 
                     data = covid)
summ(covid_online1)
covid_mademask1 <- glm(covid_mademask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                     weights = weight,
                     family = binomial(link = "probit"), 
                     data = covid)
summ(covid_mademask1)
covid_donatemoney1 <- glm(covid_donatemoney ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                       weights = weight,
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_donatemoney1)
covid_donateblood1 <- glm(covid_donateblood ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                       weights = weight,
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_donateblood1)
covid_donatetime1 <- glm(covid_donatetime ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                                                  weights = weight,
                                                  family = binomial(link = "probit"), 
                                                  data = covid)
summ(covid_donatetime1)
covid_woremask1 <- glm(covid_woremask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                          weights = weight,
                          family = binomial(link = "probit"), 
                          data = covid)
summ(covid_woremask1)
covid_volunteer1 <- glm(covid_volunteer ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, 
                       weights = weight,
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_volunteer1)

#Random effects
covid_neighbors2 <- glmer(covid_neighbors ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                        family = binomial(link = "probit"), 
                        data = covid)
summ(covid_neighbors2)
covid_neighbors2_pseudor2 <- 0.17

covid_friends2 <- glmer(covid_friends ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                      family = binomial(link = "probit"), 
                      data = covid)
summ(covid_friends2)
covid_friends2_pseudor2 <- 0.13

covid_family2 <- glmer(covid_family ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                     family = binomial(link = "probit"), 
                     data = covid)
summ(covid_family2)
covid_family2_pseudor2 <- 0.14
  
covid_online2 <- glmer(covid_online ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                     family = binomial(link = "probit"), 
                     data = covid)
summ(covid_online2)
covid_online2_pseudor2 <- 0.07

covid_mademask2 <- glmer(covid_mademask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_mademask2)
covid_mademask2_pseudor2 <- 0.08

covid_donatemoney2 <- glmer(covid_donatemoney ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                          family = binomial(link = "probit"), 
                          data = covid)
summ(covid_donatemoney2)
covid_donatemoney2_pseudor2 <- 0.15

covid_donateblood2 <- glmer(covid_donateblood ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                          family = binomial(link = "probit"), 
                          data = covid)
summ(covid_donateblood2)
covid_donateblood2_pseudor2 <- 0.14
  
covid_donatetime2 <- glmer(covid_donatetime ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                         family = binomial(link = "probit"), 
                         data = covid)
summ(covid_donatetime2)
covid_donatetime2_pseudor2 <- 0.07

covid_woremask2 <- glmer(covid_woremask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_woremask2)
covid_woremask2_pseudor2 <- 0.28

covid_volunteer2 <- glmer(covid_volunteer ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), 
                        family = binomial(link = "probit"), 
                        data = covid)
summ(covid_volunteer2)
covid_volunteer2_pseudor2 <- 0.12

##Table A14: Pandemic Actions (Probit)
stargazer(covid_neighbors1, covid_friends1, covid_family1, covid_online1, covid_mademask1, covid_donatemoney1, covid_donateblood1, covid_donatetime1, covid_woremask1, covid_volunteer1, 
          type="latex", out="Tables/covid_probit.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Neighbors}", "\\thead{Friends}", "\\thead{Family}", "\\thead{Online Groups}", "\\thead{Made Mask}", "\\thead{Donate Money}", "\\thead{Donate Blood}", "\\thead{Donate Time}", "\\thead{Wore Mask}", "\\thead{In-Person}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"))
)

##Table A15: Pandemic Actions (Random Effects)
stargazer(covid_neighbors2, covid_friends2, covid_family2, covid_online2, covid_mademask2, covid_donatemoney2, covid_donateblood2, covid_donatetime2, covid_woremask2, covid_volunteer2, 
          type="latex", out="Tables/covid_random.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Neighbors}", "\\thead{Friends}", "\\thead{Family}", "\\thead{Online Groups}", "\\thead{Made Mask}", "\\thead{Donate Money}", "\\thead{Donate Blood}", "\\thead{Donate Time}", "\\thead{Wore Mask}", "\\thead{In-Person}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                         c("Pseudo R$^{2}$", covid_neighbors2_pseudor2, covid_friends2_pseudor2, covid_family2_pseudor2, covid_online2_pseudor2, covid_mademask2_pseudor2, covid_donatemoney2_pseudor2, covid_donateblood2_pseudor2, covid_donatetime2_pseudor2, covid_woremask2_pseudor2, covid_volunteer2_pseudor2))
)

##Figure -- Marginal effects for pandemic actions
##Produce marginal effects for key IVs for each model
marg_covid_woremask1 <- summary(marginaleffects(covid_woremask1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_covid_woremask1$type <- "covid_woremask1"
marg_covid_neighbors1 <- summary(marginaleffects(covid_neighbors1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_covid_neighbors1$type <- "covid_neighbors1"
marg_covid_family1 <- summary(marginaleffects(covid_family1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_covid_family1$type <- "covid_family1"
marg_covid_donatetime1 <- summary(marginaleffects(covid_donatetime1)) %>%
  filter(term=="indiv_moral_irt_group" | term=="indiv_econ" | term=="indiv_horiz" | term=="indiv_vert" | term=="ideo5_new" | (term=="pid7" & contrast=="Strong Republican - Independent"))
marg_covid_donatetime1$type <- "covid_donatetime1"

marg_covid <- rbind(marg_covid_woremask1, marg_covid_neighbors1, marg_covid_family1, marg_covid_donatetime1)
rm(marg_covid_woremask1, marg_covid_neighbors1, marg_covid_family1, marg_covid_donatetime1)
marg_covid <- marg_covid %>% mutate(measure = case_when(term=="indiv_moral_irt_group" ~ "Moral",
                                                          term=="indiv_econ" ~ "Economic",
                                                          term=="indiv_horiz" ~ "Horizontal",
                                                          term=="indiv_vert" ~ "Vertical",
                                                          term=="ideo5_new" ~ "Ideology",
                                                          term=="pid7" ~ "Partisanship"))
marg_covid <- marg_covid %>% mutate(dv = case_when(type=="covid_woremask1" ~ "Wore Mask",
                                                     type=="covid_neighbors1" ~ "Help Neighbors",
                                                     type=="covid_family1" ~ "Help Family",
                                                     type=="covid_donatetime1" ~ "Give Time"))
marg_covid$measure <- factor(marg_covid$measure, levels=c("Moral", "Economic", "Horizontal", "Vertical", "Ideology", "Partisanship"))

##Figure 6: Effects of Individualism on Willingness to Take Collective Action
fig_pan_actions <- ggplot(data = marg_covid, aes(x=dv, y=estimate, ymin=estimate-std.error*qnorm(.975), ymax=estimate+std.error*qnorm(.975), shape=as.factor(measure))) + 
  expand_limits(y=c(-.50,.50)) + 
  geom_hline(aes(yintercept=0), linetype="dashed", color="red")+
  geom_point(size=4, stat="identity", position=position_dodge(0.4)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(0.4)) +
  ylab("Change in Probability of Support") +
  xlab("Dependent Variable") +
  coord_flip()+
  facet_wrap(~measure)+
  theme_bw()+
  theme(legend.position = "none")
fig_pan_actions
#ggsave("Figures/fig_forced_choice.pdf", width = 10.57, height = 9.55, dpi = 320)


##Summary tables of collective actions
covid_numact1 <- lm(covid_numact ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(covid_numact1)
covid_numact2 <- lmer(covid_numact ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1 | referent), data=covid)
summ(covid_numact2)
covid_numact2_pseudor2 <- 0.16
pol_numact1 <- lm(pol_act ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(pol_numact1)
pol_numact2 <- lmer(pol_act ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1 | referent), data=covid)
summ(pol_numact2)
pol_numact2_pseudor2 <- 0.17
vol_numact1 <- lm(vol_act ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=covid, weights = weight)
summ(vol_numact1)
vol_numact2 <- lmer(vol_act ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=covid)
summ(vol_numact2)
vol_numact2_pseudor2 <- 0.09

##Details for computing random effects models with different packages, discused in Footnote 15.
##Alternative method for random effects
pol_numact2_alt <- plm(pol_act ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + pid7 + age + female + race_nonewhite + college_grad + faminc_quintile, 
                        index="referent",
                        model="random",
                        effect="twoway",
                        random.method="walhus",
                        #random.models = "pooling",
                        weights = weight,
                        data=covid)
summary(pol_numact2_alt)


##RStata for random effects models

stata("random_effects_analysis.do")
test_rstata <- as.list(read.dta13("pol_act.dta"))
stargazer(test_rstata,
          type="latex")

##Table 6 Collective Action (OLS/Random Effects), 2020 YouGov Survey
stargazer(covid_numact1, covid_numact2, pol_numact1, pol_numact2, vol_numact1, vol_numact2, 
          type="latex", out="Tables/action_summary_reduced.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Covid Assistance}", "\\thead{Political Activity}", "\\thead{Volunteering}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", covid_numact2_pseudor2, "", pol_numact2_pseudor2, "", vol_numact2_pseudor2))
)

##Table A16: Collective Action (OLS/Random Effects), Full Results with Demographic Controls
stargazer(covid_numact1, covid_numact2, pol_numact1, pol_numact2, vol_numact1, vol_numact2, 
          type="latex", out="Tables/action_summary_full.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Covid Assistance}", "\\thead{Political Activity}", "\\thead{Volunteering}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", covid_numact2_pseudor2, "", pol_numact2_pseudor2, "", vol_numact2_pseudor2))
)


##Robustness checks
#Add controls for confidence in Trump
covid$trump_conf <- (4-as.numeric(covid$q46_6))/3


covid_numact3 <- lm(covid_numact ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf, data=covid, weights = weight)
summ(covid_numact3)
covid_numact4 <- lmer(covid_numact ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf + (1|referent), data=covid)
summ(covid_numact4)
covid_numact4_pseudor2 <- 0.17
covid_woremask3 <- glm(covid_woremask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf, 
                       weights = weight,
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_woremask3)
covid_woremask3_pseudor2 <- 0.15
covid_woremask4 <- glmer(covid_woremask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf + (1|referent), 
                         family = binomial(link = "probit"), 
                         data = covid)
summ(covid_woremask4)
covid_woremask4_pseudor2 <- 0.30

#Add controls for prosociality!
covid_numact5 <- lm(covid_numact ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf + prosocial, data=covid, weights = weight)
summ(covid_numact5)
covid_numact6 <- lmer(covid_numact ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf + prosocial + (1|referent), data=covid)
summ(covid_numact6)
covid_numact6_pseudor2 <- 0.23
covid_woremask5 <- glm(covid_woremask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf + prosocial, 
                       weights = weight,
                       family = binomial(link = "probit"), 
                       data = covid)
summ(covid_woremask5)
covid_woremask5_pseudor2 <- 0.17
covid_woremask6 <- glmer(covid_woremask ~ indiv_moral_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + trump_conf + prosocial + (1|referent), 
                         family = binomial(link = "probit"), 
                         data = covid)
summ(covid_woremask6)
covid_woremask6_pseudor2 <- 0.33  

##Table A17: Pandemic Collective Action (OLS/Random Effects), Controls for Trump Support and Prosociality
stargazer(covid_numact3, covid_numact4, covid_numact5, covid_numact6, covid_woremask3, covid_woremask4, covid_woremask5, covid_woremask6, 
          type="latex", out="Tables/action_summary_full_robust.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income", "Confidence in Trump", "Prosociality"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Covid Assistance}", "\\thead{Wore Mask}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", covid_numact4_pseudor2, "", covid_numact6_pseudor2, covid_woremask3_pseudor2, covid_woremask4_pseudor2, covid_woremask5_pseudor2, covid_woremask6_pseudor2))
)


##Table A18: Collective Action Pre-Covid, 2018 YouGov Survey
##2018 Volunteering and Political Participation
yougov_polact1 <- lm(part_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights = weight)
summ(yougov_polact1)
yougov_polact2 <- lmer(part_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile + (1|CSED100_2), data=yougov, weights = weight)
summ(yougov_polact2)
yougov_polact2_pseudor2 <- 0.19
yougov_volact1 <- lm(vol_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights = weight)
summ(yougov_volact1)
yougov_volact2 <- lmer(vol_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile + (1|CSED100_2), data=yougov, weights = weight)
summary(yougov_volact2)
summ(yougov_volact2)
yougov_volact2_pseudor2 <- 0.09

stargazer(yougov_polact1, yougov_polact2, yougov_volact1, yougov_volact2,
          type="latex", out="Tables/yougov_collective.tex", float = FALSE, #style="apsr",
          #column.labels   = c("\\thead{Liberty vs. Security}", "\\thead{Freedom vs. Death}"),
          #column.separate = c(2, 2), 
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Political Activity}", "\\thead{Volunteering}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", yougov_polact2_pseudor2, "", yougov_volact2_pseudor2))
)


##2020 Western States Analysis
##This limits the political participation variable to the same 5 items on the YouGov and Covid surveys. Use pol_act for an expanded list, which shows similar dynamics.
western_polact1 <- lm(pol_act2 ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=western, weights = weight)
summ(western_polact1)
western_polact2 <- lmer(pol_act2 ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=western, weights = weight)
summ(western_polact2)
western_polact2_pseudor2 <- 0.17

western_covidact1 <- lm(covid_numact ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=western, weights = weight)
summ(western_covidact1)
western_covidact2 <- lmer(covid_numact ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=western, weights = weight)
summ(western_covidact2)
western_covidact2_pseudor2 <- 0.13

#Expanded political activity index
western_polact1_full <- lm(pol_act ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=western, weights = weight)
summ(western_polact1_full)
western_polact2_full <- lmer(pol_act ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=western, weights = weight)
summ(western_polact2_full)
western_polact2_full_pseudor2 <- 0.19


##Table A19: Collective Action, 2020 Western States Survey
##Table without the vaccine results
stargazer(western_covidact1, western_covidact2, western_polact1, western_polact2, western_polact1_full, western_polact2_full, 
          type="latex", out="Tables/western_collective.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Covid Assistance}", "\\thead{Political Activity}", "\\thead{Expanded \\ Political Activity}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                         c("Pseudo R$^{2}$", "", western_covidact2_pseudor2, "", western_polact2_pseudor2, "", western_polact2_full_pseudor2))
)


##Probability of Getting Vaccinated
western_vaccine1 <- lm(vaccine ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile, data=western, weights = weight)
summ(western_vaccine1)
western_vaccine2 <- lmer(vaccine ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=western)
summ(western_vaccine2)
western_vaccine2_pseudor2 <- 0.11
western_vaccine3 <- lmer(vaccine ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + as.factor(pid7) + age + female + race_nonewhite + college_grad + faminc_quintile + (1|referent), data=western, weights = weight)
summ(western_vaccine3)
western_vaccine3_pseudor2 <- 0.12

##Alternative method for random effects
western_vaccine2_alt <- plm(vaccine ~ indiv_moral_irt_group + indiv_econ + authoritar + ideo5_new + pid7 + age + female + race_nonewhite + college_grad + faminc_quintile, 
                       index="referent",
                       model="random",
                       effect="individual",
                       random.method="walhus",
                       #random.models = "pooling",
                       weights = weight,
                       data=western)
summary(western_vaccine2_alt)

western_stata_vax <- 'use "/Users/cfk/Dropbox/Individualism Project/Data/Western States Survey/WesternStates2020_working.dta", clear
xtset referent
xtreg vaccine indiv_moral_irt_group indiv_econ authoritar ideo5_new b4.pid7 age female race_nonewhite college_grad faminc2, re
'
test_western_rstata <- stata(western_stata_vax)

##Table A20: Probability of Getting Vaccinated, 2020 Western States Survey
##Western Vaccine Table (Stata results added later by hand)
stargazer(western_vaccine1, western_vaccine3, western_vaccine2_alt, 
          type="latex", out="Tables/western_vaccine.tex", float = FALSE, #style="apsr",
          column.labels   = c("\\thead{OLS}", "\\thead{Random Effects\\\ \\texttt{LME4}}", "\\thead{Random Effects \\\ \\texttt{PLM}}"),
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Authoritarianism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Self-Reported Probability of Getting Vaccinated}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "Yes", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", western_vaccine3_pseudor2, "", ""))
)


##2021 and 2022 UNC vaccine analysis
unc_vaccine_w6_1 <- glm(vax_w6 ~ indiv_moral_irt_group_w6 + ideo5_new_w6 + pid7_w6 + as.numeric(age_w6) + female_w6 + race_nonwhite_w6 + college_grad_w6 + faminc_quintile_w6, family = binomial(link = "probit"), data=unc)
summ(unc_vaccine_w6_1)
unc_vaccine_w6_1_pseudor2 <- 0.14
unc_vaccine_w6_2 <- glmer(vax_w6 ~ indiv_moral_irt_group_w6 + ideo5_new_w6 + pid7_w6 + as.numeric(age_w6) + female_w6 + race_nonwhite_w6 + college_grad_w6 + faminc_quintile_w6 + (1|referent_w6), family = binomial(link = "probit"), data=unc)
summ(unc_vaccine_w6_2)
unc_vaccine_w6_2_pseudor2 <- 0.27

unc_vaccine_w7_1 <- glm(vax_w7 ~ indiv_moral_irt_group_w7 + ideo5_new_w7 + pid7_w7 + as.numeric(age_w7) + female_w7 + race_nonwhite_w7 + college_grad_w7 + faminc_quintile_w7, family = binomial(link = "probit"), data=unc)
summ(unc_vaccine_w7_1)
unc_vaccine_w7_1_pseudor2 <- 0.16
unc_vaccine_w7_2 <- glmer(vax_w7 ~ indiv_moral_irt_group_w7 + ideo5_new_w7 + pid7_w7 + as.numeric(age_w7) + female_w7 + race_nonwhite_w7 + college_grad_w7 + faminc_quintile_w7 + (1|referent_w6), family = binomial(link = "probit"), data=unc)
summ(unc_vaccine_w7_2)
unc_vaccine_w7_2_pseudor2 <- 0.34


##Figure 7: Effects of Individualism on Self-Reported Vaccination Status
marg_unc_vaccine_w6 <- summary(marginaleffects(unc_vaccine_w6_1)) %>%
  filter(term=="indiv_moral_irt_group_w6" | term=="ideo5_new_w6" | (term=="pid7_w6" & contrast=="Strong Republican - Independent"))
marg_unc_vaccine_w6$type <- "2021"
marg_unc_vaccine_w7 <- summary(marginaleffects(unc_vaccine_w7_1)) %>%
  filter(term=="indiv_moral_irt_group_w7" | term=="ideo5_new_w7" | (term=="pid7_w7" & contrast=="Strong Republican - Independent"))
marg_unc_vaccine_w7$type <- "2022"

marg_unc_vax <- rbind(marg_unc_vaccine_w6, marg_unc_vaccine_w7)
rm(marg_unc_vaccine_w6, marg_unc_vaccine_w7)
marg_unc_vax <- marg_unc_vax %>% mutate(measure = case_when(term=="indiv_moral_irt_group_w6" ~ "Moral Individualism",
                                                            term=="indiv_moral_irt_group_w7" ~ "Moral Individualism",
                                                            term=="ideo5_new_w6" ~ "Ideology",
                                                            term=="ideo5_new_w7" ~ "Ideology",
                                                            term=="pid7_w6" ~ "Partisanship",
                                                            term=="pid7_w7" ~ "Partisanship"))
marg_unc_vax <- marg_unc_vax %>% mutate(dv = case_when(type=="2021" ~ "2021 Respondents",
                                                   type=="2022" ~ "2022 Respondents"))
marg_unc_vax$measure <- factor(marg_unc_vax$measure, levels=c("Moral Individualism", "Ideology", "Partisanship"))

fig_vax <- ggplot(data = marg_unc_vax, aes(x=dv, y=estimate, ymin=estimate-std.error*qnorm(.975), ymax=estimate+std.error*qnorm(.975), shape=as.factor(measure))) + 
  expand_limits(y=c(-.50,.50)) + 
  geom_hline(aes(yintercept=0), linetype="dashed", color="red")+
  geom_point(size=4, stat="identity", position=position_dodge(0.4)) +
  geom_errorbar(aes(), width=.1, position=position_dodge(0.4)) +
  ylab("Change in Probability of Reporting Vaccination") +
  xlab("Survey Wave") +
  coord_flip()+
  scale_shape_manual(values=c(16, 7, 8))+
  #facet_wrap(~type)+
  theme_bw()+
  theme(legend.position = "bottom", legend.title=element_blank())
fig_vax
#ggsave("Figures/fig_forced_choice.pdf", width = 10.57, height = 9.55, dpi = 320)

##Table A21: Self-Reported Vaccination Status, 2021 and 2022 UNC Survey
#Supporting Table
#Reformatted by hand
stargazer(unc_vaccine_w6_1, unc_vaccine_w6_2, unc_vaccine_w7_1, unc_vaccine_w7_2,  
          type="latex", out="Tables/unc_vaccine.tex", float = FALSE, #style="apsr",
          covariate.labels = c("Moral Individualism", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{2021 Vaccine Status}", "\\thead{2022 Vaccine Status}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "No", "No", "No", "No"),
                         c("Pseudo R$^{2}$", unc_vaccine_w6_1_pseudor2, unc_vaccine_w6_2_pseudor2, unc_vaccine_w7_1_pseudor2, unc_vaccine_w7_2_pseudor2))
)


##Relationship to Schwartz indicators

##Figure A6: Correlations between Moral Individualism and Proxy Measures of Schwartz Values, 2018 YouGov Survey
##Correlation matrix
corrs_schwartz <- c("indiv_self_irt_group", "self_direct", "conformity", "tradition", "benevolence", "universalism")

corrs_schwartz_yougov <- yougov[corrs_schwartz]
lowerCor(corrs_schwartz_yougov)
#plot(threefactor)
pairs.panels(corrs_schwartz_yougov,pch='.', cex.labels=0.6, labels=c("Moral", "Self-Direction", "Conformity", "Tradition", "Benevolence", "Universalism"))

##Table A22: 2018 Political Activity and Volunteering with Controls for Proxy Measures of Schwartz Values
yougov_polact_schwartz <- lm(part_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + self_direct + conformity + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights = weight)
summ(yougov_polact_schwartz)
yougov_polact2_schwartz <- lmer(part_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + self_direct + conformity + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile + (1|CSED100_2), data=yougov, weights = weight)
summ(yougov_polact2)
yougov_polact2_schwartz_pseudor2 <- 0.20
yougov_volact_schwartz <- lm(vol_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + self_direct + conformity + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile, data=yougov, weights = weight)
summ(yougov_volact_schwartz)
yougov_volact2_schwartz <- lmer(vol_index ~ indiv_self_irt_group + indiv_econ + indiv_horiz + indiv_vert + authoritar + self_direct + conformity + ideo5_new + as.factor(pid7) + age + female + race_nonwhite + college_grad + faminc_quintile + (1|CSED100_2), data=yougov, weights = weight)
summ(yougov_volact2_schwartz)
yougov_volact2_schwartz_pseudor2 <- 0.10

stargazer(yougov_polact_schwartz, yougov_polact2_schwartz, yougov_volact_schwartz, yougov_volact2_schwartz,  
          type="latex", out="Tables/schwartz.tex", float = FALSE, #style="apsr",
          column.labels   = c("\\thead{OLS}", "\\thead{Random Effects}", "\\thead{OLS}", "\\thead{Random Effects}"),
          #column.separate = c(2, 2), 
          #omit = c("age", "female", "race_nonewhite", "college_grad", "faminc_quintile"),
          covariate.labels = c("Moral Individualism", "Economic Individualism", "Horizontal Individualism", "Vertical Individualism", "Authoritarianism", "Self-Direction", "Conformity", "Ideology", "Strong Democrat", "Not very strong Democrat", "Lean Democrat", "Lean Republican", "Not very strong Republican", "Strong Republican", "Age", "Female", "Nonwhite", "College Graduate", "Family Income"),
          omit.stat=c("f", "ser", "rsq", "ll", "aic", "bic"), dep.var.labels=c("\\thead{Political Activity}", "\\thead{Volunteering}"), 
          omit.table.layout = "n",
          no.space = TRUE,
          digits = 2,
          add.lines=list(c("Survey Weights", "Yes", "No", "Yes", "No"),
                         c("Pseudo R$^{2}$", "", yougov_polact2_schwartz_pseudor2, "", yougov_volact2_schwartz_pseudor2))
)

